% Smets-Wouters Model with Blanchard-Yaari Households
% augmented with anticipated monetary shocks
%
% 10/11/2015
% based on CN_model_CP_MG.m in: Dropbox/DSGE/Long_Rate/Nistico_model/CN_Model

clear all; close all;
%% ========================================================================
%%  choose model
%% ========================================================================
mdl = 8032;%9072;
wksp = 0 %41; %0
% 0: original; 2: p=0; 3: p=3; 4: estimated 
% Load model 9971's parameters

gpath = ['p0 est Fx beta baseline ',num2str(mdl),' wkspace ',num2str(wksp)];
gpath_show = []; %% gpath

if ~exist(gpath,'dir'); mkdir(gpath); end;


 load(['workspace_8031_SW.mat'])

outparatxt = [gpath '/Params.txt'];
fid01 = fopen(outparatxt,'w');



% 2015-10-09 MJS: What I did to generate workspace_907
% spec_907_704_2;
% infile0='/data/dsge_data_dir/dsgedata/save/Euro_misdata/mhparam90723547040NaN20152520Mode';
% load_modal_params;
% save('workspace_907.mat','para','mspec','subspec','dataset','stime','dates','infile0','zerobound','nant')

nant = 8; %8 % number of anticipated monetary shocks

eval(['states',num2str(mdl)]);
eval(['eqs',num2str(mdl)]);
nstates = n_end+n_exo+n_exp;
neq = nstates;

nshocks = nex-nant;
nsim = 16%2*nant   % number of periods for simulation (irfs)

% Equilibrium conditions:
% G0 z(t) = GC + G1 z(t-1) + PSI e(t) + PIE eta(t)
G0 = zeros(neq,nstates);
G1 = zeros(neq,nstates);
GC =  zeros(neq,1);
PSI = zeros(neq,nex);
PIE =  zeros(neq,nend);


%% ========================================================================
% Loop on p (death probability)
%  ========================================================================

% p = 0   : model with infinitely lived households
% p = 0.03: Nistico (2012) calibration
% p = 0.1573: model calibration implied by McKay Nakamura and Steinsson
p_ser = [0, .03, .06]%, .1292]; %
if exist('p_est')
  p_ser = [p_ser p_est];
end

zt = zeros(nstates,nsim,nant+1,size(p_ser,2));
ztpeg = zeros(nstates,nsim,nant+1,size(p_ser,2));
ztpegone = zeros(nstates,nsim,nant+1,size(p_ser,2));
ztpegoner = zeros(nstates,nsim,nant+1,size(p_ser,2));

ip=1;
for p = p_ser
  
  if ( mdl == 9072 ) 

    % Additional parameters
    para(23) = p;
    
    % Compute steady state with new p
    [alp,zeta_p,iota_p,del,ups,Bigphi,s2,h,ppsi,nu_l,zeta_w,iota_w,law,laf,betil,Rstarn,psi1,psi2,psi3,pistar,sigmac,rho,epsp,epsw...
     gam,Lmean,Lstar,gstar,rho_g,rho_b,rho_mu,rho_z,rho_laf,rho_law,rho_rm,rho_sigw,rho_mue,rho_gamm,rho_pist...
     sig_g,sig_b,sig_mu,sig_z,sig_laf,sig_law,sig_rm,sig_sigw,sig_mue,sig_gamm,sig_pist,eta_gz,eta_laf,eta_law...
     zstar,rstar,rkstar,wstar,nu_log,cstar,kstar,kbarstar,istar,ystar,sprd,zeta_spb,gammstar,vstar,nstar,...
     zeta_nRk,zeta_nR,zeta_nsigw,zeta_spsigw,zeta_nmue,zeta_spmue,zeta_nqk,zeta_nn,p,s_c,tstar,betag,etabet,bgstar] = getpara00_9072(para);
    
  elseif ( mdl == 8032 ) 
para(19) = p;

    [alp,zeta_p,iota_p,del,ups,Bigphi,s2,h,ppsi,nu_l,zeta_w,iota_w,law,laf,betil,Rstarn,psi1,psi2,psi3,pistar,sigmac,rho,epsp,epsw...
    gam,Lmean,Lstar,gstar,rho_g,rho_b,rho_mu,rho_z,rho_laf,rho_law,rho_rm,...
    sig_g,sig_b,sig_mu,sig_z,sig_laf,sig_law,sig_rm,eta_gz,eta_laf,eta_law...
    zstar,rstar,rkstar,wstar,nu_log,cstar,kstar,kbarstar,istar,ystar,p,s_c,tstar,betag,etabet,bgstar] = getpara00_8032(para);
  
  elseif ( mdl == 8031 ) 
para(19) = p;
  

    [alp,zeta_p,iota_p,del,ups,Bigphi,s2,h,ppsi,nu_l,zeta_w,iota_w,law,laf,betil,Rstarn,psi1,psi2,psi3,pistar,sigmac,rho,epsp,epsw...
    gam,Lmean,Lstar,gstar,rho_g,rho_b,rho_mu,rho_z,rho_laf,rho_law,rho_rm,...
    sig_g,sig_b,sig_mu,sig_z,sig_laf,sig_law,sig_rm,eta_gz,eta_laf,eta_law...
    zstar,rstar,rkstar,wstar,nu_log,cstar,kstar,kbarstar,istar,ystar,p,s_c,tstar,betag,etabet,bgstar] = getpara00_8031(para);
  
  elseif ( mdl == 803 ) 
    para([19 20]) = [];
    
    [alp,zeta_p,iota_p,del,ups,Bigphi,s2,h,ppsi,nu_l,zeta_w,iota_w,law,laf,bet,Rstarn,psi1,psi2,psi3,pistar,sigmac,rho,epsp,epsw...
    gam,Lmean,Lstar,gstar,rho_g,rho_b,rho_mu,rho_z,rho_laf,rho_law,rho_rm,...
    sig_g,sig_b,sig_mu,sig_z,sig_laf,sig_law,sig_rm,eta_gz,eta_laf,eta_law...
    zstar,rstar,rkstar,wstar,wl_c,cstar,kstar,kbarstar,istar,ystar] = getpara00_803(para);
    
  end
  
  
  %% write params
  fprintf(fid01,'\n\n p = %3.4f  \n',[p]);
  if exist('etabet');   
    fprintf(fid01,' eta = %3.4f; ',[etabet]); 
  end;
  if exist('betil');   
    kappa = ((1-zeta_p*betil*exp((1-sigmac)*zstar))*(1-zeta_p))/(zeta_p*((Bigphi-1)*epsp+1))*1/(1+iota_p*betil*exp((1-sigmac)*zstar));

    fprintf(fid01,' betil/(1+betil) = %3.4f; ',[betil/(1+betil)]); 
    fprintf(fid01,' kappa  = %3.4f; ',[kappa]); 
    fprintf(fid01,' betil/(1+iota_p*betil) = %3.4f; ',[betil/(1+iota_p*betil)]); 
  else
    kappa = ((1-zeta_p*bet*exp((1-sigmac)*zstar))*(1-zeta_p))/(zeta_p*((Bigphi-1)*epsp+1))*1/(1+iota_p*bet*exp((1-sigmac)*zstar));

    fprintf(fid01,' bet/(1+bet) = %3.4f; ',[bet/(1+bet)]); 
    fprintf(fid01,' kappa  = %3.4f; ',[kappa]); 
    fprintf(fid01,' bet/(1+iota_p*bet) = %3.4f; ',[bet/(1+iota_p*bet)]); 
    
  end;
  
  
  % Build the model's matrices
  eval(['eqcond',num2str(mdl)]);

  
  % Steady state
  zss = zeros(nstates,1);
  %zss = (G0-G1)\GC;

  % Solving using gensys. Solution is of the form:
  % z(t) = Ch+Th*z(t-1)+Rh*e(t)+ywt*inv(I-fmat*inv(L))*fwt*e(t+1) 
  [Th,Ch,Rh,fmat,fwt,ywt,gev,eu] = gensys(G0,G1,GC,PSI,PIE);  
  

  %% ========================================================================
  %%   Simulate the model

    
  for H = 0:nant

    % IRFs to anticipated shock at horizon H
    % --------------------------------------
    et = zeros(nex,nsim);
    %et(iepsi,1) = -0.0025/4; % unant. pol shock: reduction of 25bps annualized
    %et(7+H,1) = -0.0025/4; %  ant. pol shock in 4Q: reduction of 25bps annualized
    if H == 0
        et(rm_sh,1)= -0.0025/4;
    else
        eval(strcat('et(rm_shl',num2str(H),',1)= -0.0025/4;')); %  ant. pol shock in H Q: reduction of 25bps annualized    
    end
    zt1 = zss;
    for t = 1:nsim
        zt(:,t,H+1,ip) =  Th*zt1 + Rh*et(:,t); %Ch +
        zt1 = zt(:,t,H+1,ip);
    end    


    % IRFs to FFR peg for H periods
    % -----------------------------
    % Set anticipated shocks to peg FFR from t+1 to t+H
    FFRpeg = -0.25/400;
    %% constant
    PsiR1 = 0;
    %% ZZ matrix
    PsiR2 = zeros(1,nstates); PsiR2(1,R_t)=1;
    if H == 0
        Rht = Rh(:,rm_sh);
    else        
        %Rht = Rh(:,[iepsi+1:iepsi+H]); % Columns of Rh referring to anticipated monetary shocks
        eval(strcat('Rht = Rh(:,[rm_sh,rm_shl1:rm_shl',num2str(H),']);')); % Columns of Rh referring to anticipated monetary shocks
    end
    stinit = zss; % initial state vector
    bb = zeros(H+1,1);
    MH = zeros(H+1,H+1);
    for hh = 1:H+1;      
      bb(hh,1) = FFRpeg - PsiR1 - PsiR2*(Th)^hh*stinit;
      MH(hh,:) = PsiR2*(Th)^(hh-1)*Rht;
    end
    monshocks = MH\bb;
    
    % Simulate the model
    etpeg = zeros(nex,nsim);
    if H == 0
        etpeg(rm_sh,1) = monshocks;
    else       
        %etpeg(iepsi,1) = monshocks(1); % unant. pol shock
        %etpeg(iepsi:iepsi+H-1,1) = monshocks; %  ant. pol shock
        eval(strcat('etpeg([rm_sh,rm_shl1:rm_shl',num2str(H),'],1) = monshocks; ')); 
    end
    ztpeg1 = zss;
    for t = 1:nsim
      ztpeg(:,t,H+1,ip) = Ch + Th*ztpeg1 + Rh*etpeg(:,t);
      ztpeg1 = ztpeg(:,t,H+1,ip);
    end
    

    % IRFs to FFR peg for 1 period, H periods ahead, and fixed to zero otherwise
    % -----------------------------
    % Set anticipated shocks to peg FFR from t+1 to t+H
    FFRpeg = -0.25/400;
    %% constant
    PsiR1 = 0;
    %% ZZ matrix
    PsiR2 = zeros(1,nstates); PsiR2(1,R_t)=1;
    if H == 0
        Rht = Rh(:,rm_sh);
    else        
        %Rht = Rh(:,[iepsi+1:iepsi+H]); % Columns of Rh referring to anticipated monetary shocks
        eval(strcat('Rht = Rh(:,[rm_sh,rm_shl1:rm_shl',num2str(H),']);')); % Columns of Rh referring to anticipated monetary shocks
    end
    stinit = zss; % initial state vector
    bb = zeros(H+1,1);
    MH = zeros(H+1,H+1);
    for hh = 1:H+1;
            if hh < H+1
        FFRpeg_ = 0;            
            else
        FFRpeg_ = FFRpeg;
            end
      bb(hh,1) = FFRpeg_ - PsiR1 - PsiR2*(Th)^hh*stinit;
      MH(hh,:) = PsiR2*(Th)^(hh-1)*Rht;
    end

    monshocks = MH\bb;
    
    % Simulate the model
    etpeg1 = zeros(nex,nsim);
    if H == 0
        etpeg1(rm_sh,1) = monshocks;
    else       
        %etpeg1(iepsi,1) = monshocks(1); % unant. pol shock
        %etpeg1(iepsi:iepsi+H-1,1) = monshocks; %  ant. pol shock
        eval(strcat('etpeg1([rm_sh,rm_shl1:rm_shl',num2str(H),'],1) = monshocks; ')); 
    end
    ztpegone1 = zss;
    for t = 1:nsim
      ztpegone(:,t,H+1,ip) = Ch + Th*ztpegone1 + Rh*etpeg1(:,t);
      ztpegone1 = ztpegone(:,t,H+1,ip);
    end
    
    % IRFs to REAL RATE peg for 1 period, H periods ahead, and fixed to zero otherwise
    % -----------------------------
    % Set anticipated shocks to peg FFR from t+1 to t+H
    FFRpeg = -0.25/400;
    %% constant
    PsiR1 = 0;
    %% ZZ matrix
    PsiR2 = zeros(1,nstates); PsiR2(1,R_t)=1; PsiR2(1,E_pi)= -1;
    if H == 0
        Rht = Rh(:,rm_sh);
    else        
        %Rht = Rh(:,[iepsi+1:iepsi+H]); % Columns of Rh referring to anticipated monetary shocks
        eval(strcat('Rht = Rh(:,[rm_sh,rm_shl1:rm_shl',num2str(H),']);')); % Columns of Rh referring to anticipated monetary shocks
    end
    stinit = zss; % initial state vector
    bb = zeros(H+1,1);
    MH = zeros(H+1,H+1);
    for hh = 1:H+1;
            if hh < H+1
        FFRpeg_ = 0;            
            else
        FFRpeg_ = FFRpeg;
            end
      bb(hh,1) = FFRpeg_ - PsiR1 - PsiR2*(Th)^hh*stinit;
      MH(hh,:) = PsiR2*(Th)^(hh-1)*Rht;
    end

    monshocks = MH\bb;
    
    % Simulate the model
    etpeg1r = zeros(nex,nsim);
    if H == 0
        etpeg1r(rm_sh,1) = monshocks;
    else       
        %etpeg(iepsi,1) = monshocks(1); % unant. pol shock
        %etpeg(iepsi:iepsi+H-1,1) = monshocks; %  ant. pol shock
        eval(strcat('etpeg1r([rm_sh,rm_shl1:rm_shl',num2str(H),'],1) = monshocks; ')); 
    end
    ztpegoner1 = zss;
    for t = 1:nsim
      ztpegoner(:,t,H+1,ip) = Ch + Th*ztpegoner1 + Rh*etpeg1r(:,t);
      ztpegoner1 = ztpegoner(:,t,H+1,ip);
    end
  
  end


  %% ========================================================================
  % Figures
  

  % Effects of anticipated shocks at horizons H = 1 to nant

  H=0;
  figure(1000+ip) 
  subplot(221), plot(0:nsim-1, 400*[zt(R_t,:,H+1,ip)])
  title(['Interest rate (p = ' num2str(p) ', \eta = ' num2str(etabet) ')'])
  % axis([0 nsim -0.1 4*zss(ii)+0.1])
  subplot(222), plot(0:nsim-1, 400*[zt(pi_t,:,H+1,ip)], 0:nsim-1, zeros(1,nsim),':k')
  title(['Inflation (% annualized)'])
  subplot(223), plot(0:nsim-1, 400*[zt(c_t,:,H+1,ip)])
  title(['Consumption (level, %)'])
  subplot(224), plot(0:nsim-1, 100*[zt(y_t,:,H+1,ip)], 0:nsim-1, zeros(1,nsim),':k')
  title(['Output (level, %)'])
  
  for H = 1:8
    subplot(221),
    hold on
    plot(0:nsim-1, 400*[zt(R_t,:,H+1,ip)]);
    subplot(222),
    hold on
    plot(0:nsim-1, 400*[zt(pi_t,:,H+1,ip)]);
    subplot(223),
    hold on
    plot(0:nsim-1, 400*[zt(c_t,:,H+1,ip)]);
    subplot(224),
    hold on
    plot(0:nsim-1, 100*[zt(y_t,:,H+1,ip)]);
  end
  hold off
  
  
  % Effects of interest rate peg until horizons H = 1 to nant

  H=1;
  figure(2000+ip)
  subplot(221), plot(0:nsim-1, 400*[ztpeg(R_t,:,H,ip)])
  title(['Interest rate PEG (p = ' num2str(p) ', \eta = ' num2str(etabet) ')'])
  % axis([0 nsim -0.1 4*zss(ii)+0.1])
  subplot(222), plot(0:nsim-1, 400*[ztpeg(pi_t,:,H,ip)], 0:nsim-1, zeros(1,nsim),':k')
  title(['Inflation (% annualized)'])
  subplot(223), plot(0:nsim-1, 400*[ztpeg(c_t,:,H,ip)])
  title(['Consumption (level, %)'])
  subplot(224), plot(0:nsim-1, 100*[ztpeg(y_t,:,H,ip)], 0:nsim-1, zeros(1,nsim),':k')
  title(['Output (level, %)'])
  for H = 2:8
    subplot(221),
    hold on
    plot(0:nsim-1, 400*[ztpeg(R_t,:,H,ip)]);
    subplot(222),
    hold on
    plot(0:nsim-1, 400*[ztpeg(pi_t,:,H,ip)]);
    subplot(223),
    hold on
    plot(0:nsim-1, 400*[ztpeg(c_t,:,H,ip)]);
    subplot(224),
    hold on
    plot(0:nsim-1, 100*[ztpeg(y_t,:,H,ip)]);
  end
  hold off

  if length(p_ser) < 5
    
    LStyle = {'-','--','-.','-'};    
    legendInfo{ip} = ['p = ' num2str(p,3) ', \eta = ' num2str(etabet,3)];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    %% IRFs
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    %% H_: what horizons to show
    H_ = [0 nant];
    
    
    for H = H_
      figure(10000+H)
      subplot(221),
      P = plot(0:nsim-1, 400*[zt(R_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Interest rate (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 400*[zt(pi_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Inflation (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(223),
      hold on
      P = plot(0:nsim-1, 100*[zt(c_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Consumption (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[zt(y_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Output (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on


      if (ip == length(p_ser) )
        
        subplot(221), 
        legend(legendInfo)         
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        
        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show]);
      end

      figure(10010+H)

      subplot(221),
      P = plot(0:nsim-1, 100*[zt(L_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['L (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 100*[zt(mc_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Marginal Costs (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(223),
      hold on
      P = plot(0:nsim-1, 400*[zt(R_t,:,H+1,ip)-zt(E_pi,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Real rate (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[zt(s_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Wealth (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      if (ip == length(p_ser) )
        
        subplot(221), 
        
        legend(legendInfo);
        
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;        
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        
        %        L = legend(legendInfo) set(L,'OuterPosition',[0 0 .5 .5])  
                
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        
        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show]);
        
      end
      
    end    

    
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
     %% H_: what horizons to show for all following exercises
    H_p = [0 4 nant];

     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
     % FFR Peg
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    
    for H = H_p
      figure(20000+H)
      subplot(221),
      P = plot(0:nsim-1, 400*[ztpeg(R_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Interest rate (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 400*[ztpeg(pi_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Inflation (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(223),
      hold on
      P = plot(0:nsim-1, 100*[ztpeg(c_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Consumption (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[ztpeg(y_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Output (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on


      if (ip == length(p_ser) )

        subplot(221), 
        legend(legendInfo) 
        
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;        
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        
        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show]);
      end

      figure(20010+H)

      subplot(221),
      P = plot(0:nsim-1, 100*[ztpeg(L_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['L (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 100*[ztpeg(mc_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Marginal Costs (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(223),
      hold on
      P = plot(0:nsim-1, 400*[ztpeg(R_t,:,H+1,ip)-ztpeg(E_pi,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Real rate (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[ztpeg(s_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Wealth (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      if (ip == length(p_ser) )
        subplot(221), 
        legend(legendInfo) 
        
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;        
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        
        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show]);
      end
      
    end    

  
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
     % FFR peg for 1 period, H periods ahead, and fixed to zero otherwise
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    
    for H = H_p
      figure(21000+H)
      subplot(221),
      P = plot(0:nsim-1, 400*[ztpegone(R_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Interest rate (% annualized)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([-.25 .1])
          yticks(-.25:.05:.1)
          yticklabels({'0.25', '--0.2', '0.15', '--0.1', '0.05', '0', '0.05', '0.1'})
          set(gca, 'ticklabelinterpreter', 'latex')
      elseif H == 4
          ylim([-.3 .05])
          yticks(-.3:.05:.05)
          yticklabels({'--0.3', '0.25', '--0.2', '0.15', '--0.1', '0.05', '0', '0.05'})
          set(gca, 'ticklabelinterpreter', 'latex')
      elseif H == 8
          ylim([-.4 1])
          yticks(-.4:.2:1)
          yticklabels({'--0.4', '--0.2', '0', '0.2', '0.4', '0.6', '0.8', '1'})
          set(gca, 'ticklabelinterpreter', 'latex')
      end
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 400*[ztpegone(pi_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Inflation (% annualized)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([0 .07])
          yticks(0:.01:.07)
      elseif H == 4
          ylim([0 .35])
          yticks(0:.05:.35)
      elseif H == 8
          ylim([0 5])
          yticks(0:1:5)
      end
      hold on

      subplot(223),
      hold on
      P = plot(0:nsim-1, 100*[ztpegone(c_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Consumption (level, %)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([0 .14])
          yticks(0:.02:.14)
      elseif H == 4
          ylim([0 .5])
          yticks(0:.1:.5)
      elseif H == 8
          ylim([0 6])
          yticks(0:1:6)
      end
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[ztpegone(y_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Output (level, %)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([0 .14])
          yticks(0:.02:.14)
      elseif H == 4
          ylim([0 .7])
          yticks(0:.1:.7)
      elseif H == 8
          ylim([0 8])
          yticks(0:2:8)
      end
      hold on


      if (ip == length(p_ser) )

        subplot(221), 
        
        plot(0:nsim-1, zeros(1,nsim),':k');
        legend(legendInfo) 
        box on;        
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        

        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show])
      end

      figure(21010+H)

      subplot(221),
      P = plot(0:nsim-1, 100*[ztpegone(L_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['L (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 100*[ztpegone(mc_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Marginal Costs (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(223),
      hold on
      P = plot(0:nsim-1, 400*[ztpegone(R_t,:,H+1,ip)-ztpegone(E_pi,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Real rate (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[ztpegone(s_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Wealth (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      if (ip == length(p_ser) )
        subplot(221), 
        legend(legendInfo) 
        
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;        
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        
        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show]);
      end
      
    end    


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
     % REAL RATE peg for 1 period, H periods ahead, and fixed to zero otherwise
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    for H = H_p
      figure(22000+H)
      subplot(221),
      P = plot(0:nsim-1, 400*[ztpegoner(R_t,:,H+1,ip)-ztpegoner(E_pi,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Real rate (% annualized)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([-.25 .1])
          yticks(-.25:.05:.1)
          yticklabels({'0.25', '--0.2', '0.15', '--0.1', '0.05', '0', '0.05', '0.1'})
          set(gca, 'ticklabelinterpreter', 'latex')
      elseif H == 4
          ylim([-.3 .05])
          yticks(-.3:.05:.05)
          yticklabels({'--0.3', '0.25', '--0.2', '0.15', '--0.1', '0.05', '0', '0.05'})
          set(gca, 'ticklabelinterpreter', 'latex')
      elseif H == 8
          ylim([-.3 .05])
          yticks(-.3:.05:.05)
          yticklabels({'--0.3', '0.25', '--0.2', '0.15', '--0.1', '0.05', '0', '0.05'})
          set(gca, 'ticklabelinterpreter', 'latex')
      end
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 400*[ztpegoner(pi_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Inflation (% annualized)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([0 .06])
          yticks(0:.01:.06)
      elseif H == 4
          ylim([0 .2])
          yticks(0:.05:.2)
      elseif H == 8
          ylim([0 .35])
          yticks(0:.05:.35)
      end
      hold on

      subplot(223),
      hold on
      P = plot(0:nsim-1, 100*[ztpegoner(c_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Consumption (level, %)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([0 .1])
          yticks(0:.02:.1)
      elseif H == 4
          ylim([0 .2])
          yticks(0:.05:.2)
      elseif H == 8
          ylim([0 .25])
          yticks(0:.05:.25)
      end
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[ztpegoner(y_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Output (level, %)'], 'FontWeight', 'normal')
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      if H == 0
          ylim([0 .12])
          yticks(0:.02:.12)
      elseif H == 4
          ylim([0 .25])
          yticks(0:.05:.25)
      elseif H == 8
          ylim([0 .35])
          yticks(0:.05:.35)
      end
      hold on


      if (ip == length(p_ser) )

        subplot(221), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        legend(legendInfo) 
        box on;        
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        
        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show]);
      end

      figure(22010+H)

      subplot(221),
      P = plot(0:nsim-1, 100*[ztpegoner(L_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['L (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      subplot(222),
      P = plot(0:nsim-1, 100*[ztpegoner(mc_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Marginal Costs (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(223),
      P = plot(0:nsim-1, 400*[ztpegoner(R_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Interest rate (% annualized)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on

      subplot(224),
      hold on
      P = plot(0:nsim-1, 100*[ztpegoner(s_t,:,H+1,ip)]);
      if  (ip == 1 )
        title(['Wealth (level, %)'])
      end
      set(P,'LineWidth',2+((length(p_ser)-ip)/length(p_ser))^.5,'LineStyle',LStyle{ip},'Color',[1 1 1]*((length(p_ser)-ip)/length(p_ser))^.35);
      hold on
      
      if (ip == length(p_ser) )
        subplot(221), 
        legend(legendInfo) 
        
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;        
        hold on;
        subplot(222), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(223), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        subplot(224), 
        plot(0:nsim-1, zeros(1,nsim),':k');
        box on;
        hold on;
        
        %[s1,s2] = 
        subtitle(['Horizon  ' num2str(H) '  ' gpath_show]);
      end
      
    end    

  end  

  close(figure(1000+ip)); 
  close(figure(2000+ip)); 
  %% End loop over values of p
  ip = ip+1;

end  % loop over values of p

for H = H_
  figure(10000+H);
  hold off;
  saveas(figure(10000+H),[gpath '/comp_IRFs_H',num2str(H),'.eps'] ,'eps');

  figure(10010+H);
  hold off;
  saveas(figure(10010+H),[gpath '/comp_IRFs_2_H',num2str(H),'.eps'] ,'eps');

end

for H = H_p
  
  figure(20000+H);
  hold off;
  saveas(figure(20000+H),[gpath '/comp_PEG_H',num2str(H),'.eps'] ,'eps');

  figure(20010+H);
  hold off;
  saveas(figure(20010+H),[gpath '/comp_PEG_2_H',num2str(H),'.eps'] ,'eps');

  
    figure(21000+H);
  hold off;
  saveas(figure(21000+H),[gpath '/comp_PEG1_H',num2str(H),'.eps'] ,'eps');

  figure(21010+H);
  hold off;
  saveas(figure(21010+H),[gpath '/comp_PEG1_2_H',num2str(H),'.eps'] ,'eps');

  
    figure(22000+H);
  hold off;
  saveas(figure(22000+H),[gpath '/comp_PEG1REAL_H',num2str(H),'.eps'] ,'eps');

  figure(22010+H);
  hold off;
  saveas(figure(22010+H),[gpath '/comp_PEG1REAL_2_H',num2str(H),'.eps'] ,'eps');

end

fclose(fid01);

% pause;
close all;

